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ABSTRACT 

A theory is described based on resonant thermal diffusion waves in the sun that explains many details of the 
paleotemperature record for the last 5.3 million years. These include the observed periodicities, the relative strengths 
of each observed cycle, and the sudden emergence in time for the 100 thousand year cycle. Other prior work 
suggesting a link between terrestrial paleoclimate and solar luminosity variations has not provided any specific 
mechanism. The particular mechanism described here has been demonstrated empirically, although not previously 
invoked in the solar context. The theory, while not without its own unresolved issues, also lacks most of the problems 
associated with Milankovitch cycle theory. 



I. PALEOTEMPERATURE DATA 

The possibility that fluctuations in solar luminosity may be responsible for changes in global temperatures has 
not been overlooked by researchers, (Lean, 1997) although most explanations for periodicities in paleotemperatures are 
believed to involve factors unrelated to solar luminosity.(Zachos,2001) Nevertheless, some researchers have suggested 
that periodic solar variability has been the cause of global temperature cycles, with periods ranging from the 11 year 
sunspot period to cycles as long as 2000 years. (Hu, 2003) Earlier reports, however, suggest no specific mechanism for 
long term solar periodicity. 

Paleotemperatures are inferred from the 18 0/ 16 ratio found in benthic sediments and in ice cores versus depth. In 
particular, it is believed that a 5 °C decrease in global temperature corresponds to a 0.1% increase in the ls O/ le O ratio 
(i.e., delta 18 0) in benthic sediments. (Lisiecki, 2005) Fig. 1 shows a compilation by Lisiecki and Raymo(Lisiecki,2005) 
who stacked 57 sets of data on benthic sediments that cover 5.3 million years (My) before present. Figs. 2 and 3 show 
the Fourier amplitude spectra. 

II. DIFFUSION WAVES 

Diffusion waves are an established phenomenon which has been studied theoretically and experimen- 
tally, (Mandelis, 2000) and even for heat transport in magnetized plasmas. (Reynolds, 2001) They can occur for any 
process, such as heat flow when a time-dependent source is present, in which case, the waves represent propagating 
temperature fluctuations, ip(r,t), and they must satisfy the diffusion equation. Assuming a sinusoidal time depen- 
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dence, and separating ip(r,t) — (f>(r)e ZU3t mto space and time dependent parts, we find that <fi(r) satisfies: 



V 2 - K 2 = Q (r) (f) 
where Qo( r ) is the space-dependent part of a source term, and where the complex wave number k satisfies 




Note that according to eq. the imaginary (attenuation) part of the complex wave number k and the real (oscilla- 
tory) part are equal, implying equality of penetration depth L = 27t/to(k), and wavelength A = 2tt / Re(n) , so higher 
frequency waves are more severely damped. For a uniform thermal diffusivity D t and a spherically symmetric source 
eq. [1] yields the solution 



sin Re(n)r 

4>{r) = (3) 

For the more general case of an isotropic medium in which n(r) is dependent only on radial distance to the origin 
the same solution applies, but k must now be replaced by its average value between the source radius ro and the r at 
the wave location, i.e., 



«(r — ?'o) = / n(r)dr. (4) 

J To 



III. THE SOLAR INTERIOR 



Solar physicists use a Standard Solar Model (SSM) to calculate many properties of the sun's interior. According 
to the SSM, the three interior regions of interest here are: (1) the core (of radius about 25% of the solar radius r Q ) 
where fusion occurs, (2) the surrounding radiative zone (extends to about O.Tr©), and (3) the surrounding convection 
zone, where heat makes its way to the surface by the much more rapid process of convective flow. The fairly sharp 
boundary layer between the radiative and convection zones is known as the "tachocline." 

Normally, one associates diffusivity with heat flow by conduction, whereas radiation is taken as a line-of-sight 
transmission process. However, inside the sun where the mean free path for photons between emission and subsequent 
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reabsorption is very short (typically several millimeters), the diffusion description is applicable. Obviously, radiation 
is the dominant heat transfer process inside the core and radiative zones, owing to the high temperatures there, and 
the T 4 dependence on radiative emission and absorption rates. The diffusivity parameter depends on a variety of 
other quantities that govern the likelihood of photons being radiated and reabsorbed. In particular, it can be shown 
that the thermal diffusivity at any radial distance as (Miesch, 2005): D t — 16a s bT 3 /(3xp 2 Cp), where a 8 b is the 
Stefan-Boltzmann constant, T is the absolute temperature, \ IS the opacity, p is the density, and Cp is the specific 
heat. Each of the parameters in the previous equation can be found as a function of radial distance using data from 
the GONG collaboration (Christensen-Dalsgaard,1996). 

Within the SSM description, the solar core is considered to be in a state of quasistatic equilibrium changing only 
on a timescale of 30 million years. However, Grandspicrre and Agoston have noted, instabilities can be generated 
in a rotating plasma in the presence of a magnetic field, and these instabilities can give rise to thermal fluctuations, 
and heat waves that represent deviations from SSM's.(Grandpierre,2005) They further show that constraints on the 
size of the core magnetic field (Friedland, 2004), (Gough, 1998) are such that fluctuations can indeed arise, as a natural 
consequence of the nonlinearity of the MHD equations. 

These thermal fluctuations in the form of hot " bubbles" tend to expand, rise, and cool on timescales Grandspierre 
and Agoston have estimated, which are a function of their spatial extent and the magnitude of their temperature 
excess. (Grandpierre, 2005) Were the bubbles to remain small and isolated, the diffusion waves emanating from these 
fluctuations would tend to cancel each other out. However, Grandspicrre and Agoston have shown that the bubbles 
tend to merge and grow in size, and the decay time of the waves can grow even longer than 10 7 years, as the 
perturbations smooth out to include ever larger regions - in effect making the solar core a region of small amplitude 
thermal oscillations, instead of the quiescent state normally assumed. (Grandpierre, 2005) 

IV. DESCRIPTION OF THE THEORY 

We postulate small random variations in the temperature of the outer solar core (between 0.21r Q and 
O.25r ),(Ferro,2OO5) which are sufficiently smoothed out, so as to be a function only of radial distance. The specific 
noise spectrum is unimportant, except that the noise frequencies should have slowly changing phases over times on a 
scale of 10 5 y. 

Each of the noise frequencies may be thought of as a source of radially symmetric diffusion waves. When a 
wave reaches the tachocline, it is reflected back to the core as an ingoing wave. However, it should be noted that 
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reflections in the case of diffusion waves are not of the usual type, where well-defined wavefronts are involved, but only 
resemble conventional wave reflections in the high frequency or long path-length limit (Mandelis,200f ). Otherwise, 
a considerable amount of "phase smearing" prevents us from describing the reflection of discrete wavefronts. At a 
boundary, such as the tachocline, the depletion or accumulation of photons (and the associated radially symmetric 
periodic temperature fluctuation), creates an inward travelling diffusion wave. This inward wave interferes with the 
next outgoing wave produced at that frequency, and a new outward wave is generated from a superposition of the 
two waves. Again, however, the interference process is not one of simply superimposing well-defined wavefronts, 
which do not exist here. Instead, if the inward and outward waves both correspond to a maximum positive or 
negative fluctuation in photon density (i.e., have have matched phases), the resultant temperature fluctuation would 
be greatest, and hence resonance would be achieved. To achieve resonance the ingoing wave must arrive back at the 
core in phase with the next source cycle at this frequency, or more generally after n cycles. Thus, the condition for 
resonance for the nth mode is: 



where n =1,2,3,..., and 2r is the round trip distance travelled by the wave from emission to absorption back at the 
core. From eq. [5] we see that the period T n of the nth resonant mode T n = 2-7r/o; n oc 1/n 2 , so that 



with the period Ti of the lowest frequency resonant mode being twice the diffusion time td for photons to undergo 
a random walk through the radiation zone before they reach the tachocline. Fortunately, td can be calculated from 
the sun's internal density and opacity (both functions of radial distance r) using data from the SSM.(Christenscn- 
Dalsgaard,2005) One the most widely-cited estimates of td from Mitalas and Sills(Mitalas,1992) is 1.7 x 10 5 y = 170fcy. 
We have repeated their calculation, and obtained 190 ky. In what follows, we shall use a value that brackets these 
two: t d = 180 ± Why. 

Using eq. [5] we find the predicted periods T n for the nth resonant mode: 




(5) 



T n = T./n 2 



(6) 



(7) 
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Some of the periods (for n — 2, 3 and 4) given by eq. [7] correspond to reported periods in the paleotemperature 
data - see markers in fig. 2. But before we can have confidence that resonant solar diffusion waves may offer an 
explanation of periodic variations in paleotemperatures, we must address several key issues, including whether it is 
possible for the waves to attain a large degree of resonant amplification. 

Damping of diffusion waves Given that a diffusion wave for the nth resonant mode has a phase of 2irn on 
returning to the source, it will have been heavily damped, and will have an intrinsic damping: 

d n = e^ 2 ™ (8) 

The large damping of diffusion waves might seem to make the possibility of their playing a role in any resonant 
amplification process in the sun remote. For example, given a damping d n for a wave after one round trip, the next 
outgoing wave has an amplitude 1 + d n , and hence the number of cycles required for the amplitude of the nth mode 
to double is given as ln2/d„, yielding doubling times in ky for mode n: < X 2 = 2n~ 2 tde 27vn \n2. For n > 1 we find 
doubling times that are longer than the age of the sun - even if diffusion waves are reflected 1 00% from the tachocline. 
Yet the doubling times are lengthened even more, because of damping due to "blurry boundaries," since fractions of 
an outgoing wave are reflected back from different parts of the tachocline, and return to different parts of the core. 
Based on the relevant layer thicknesses, we find a phase smearing of the nth mode: Atfi ks ±27rn/10. This additional 
damping, which increases with n, makes it impossible for modes n > 4 to survive. 

Gain due to fusion amplification. There is a compensating amplification effect that more than offsets the 
enormous damping, and allows resonant diffusion waves to emerge with large amplitude in a relatively few cycles. 
According to the SSM, the rate of nuclear fusion in the core of the sun (and hence solar luminosity L) is proportional to 
T 4 - 2 for small deviations from the ambient core temperature. (Bahcall, 2001) To see the effect of this large exponent in 
amplifying diffusion waves, consider an outgoing wave for mode n initially created due to a positive core temperature 
fluctuation <^ n (ro) at radius ro = 0.21r Q . In advancing a small distance dr outward, the wave amplitude is reduced 
by dcj) = ndr. However, if the wave is absorbed by the core, it elevates the core temperature at that r by the same 
amount, and hence raises the luminosity there by a 4.2-fold greater fraction, due to the L oc T relationship. Thus, 
there is a net gain in wave amplitude over this distance dr of 3.2ndr. The amplification factor multiplies the size of 
the temperature fluctuation, and it occurs virtually instantaneously on a diffusion time scale, so that the amplification 
looks like a coherent radiation wave superimposed on a much slower moving thermal wave of the same frequency. 
As the wave proceeds outward through the outer core, it continues to produce an ever increasing amount of thermal 
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radiation, owing to the continual increase in luminosity. Upon reaching the edge of the core (0.25r Q ) the nth resonant 
mode (and the core temperature) has been amplified by a gain: 

9n = e 3 - 2ftAr (9) 

where R can be found using eqs. [2j SI EJ with ro = 0.21r Q and r = 0.25r Q . The same mechanism operates for the 
negative portion of the wave cycle, and it amplifies reductions in core temperature. This amplification mechanism is 
similar to the exponential increase in intensity in the case of lasing, except that here (1) fusion makes optical pumping 
unecessary, and (2) the spherical symmetry assures that there is no loss of waves "out the sides," as occurs in a laser. 
Evaluating eq. [9] numerically using SSM data yields g n = e 9 ' 6 ". 

Such a large gain (more than offsetting the e~ 27rn damping) cannot be taken as a realistic indicator of how large 
the wave can grow before it leaves the core, but it does suggest that a wave might be able to resonantly attain its 
maximum value after only a few cycles. It should be noted that a very similar resonance phenomenon to that described 
here - but without any fusion gain - has been seen empirically in thermal-wave cavity environments. (Shen, 1995) 

Obviously, the theory provided here, which ignores nonthermal magnetic interactions, and assumes an initial tem- 
perature fluctuation that is only a function of radial distance is highly simplified. Far more numerous non-spherically 
symmetric temperature fluctuations would not result in a similar amplification process, and hence they will become 
relatively less intense over time. However, both the nonspherically symmetric fluctuations and the magnetic interac- 
tions must be considered for a full understanding of losses and gains, so as to permit a determination of the maximum 
steady state radial wave amplitude. 

V. COMPARISON OF BETWEEN THEORY AND PALEOCLIMATE DATA 

Periods of each mode. As fig. [^markers) and table U show, the Fourier amplitude spectrum for the data shows 
peaks at the theory's predicted modes n = 2,3 and 4. The overall data with a fifth order background subtraction to 
eliminate low frequency leakage shows no indication of the n = 1 mode in fig. [2J With no background subtraction 
there is a non-statistically significant suggestion of a peak corresponding to this mode in fig.[2Ic) and especially (d), 
which together correspond to the interval 2048 to 4095 years before present. Were it present, a n = 1 peak (located 
at f = 0.28 cycles per 100 ky) could not easily be distinguished from a steeply rising low frequency background. 

No periodicities predicted by the model are clearly absent from the data. It is true, however, that one small peak 
seen in the overall spectrum (but not in the quartile spectra) with a period of 19 ky (dotted arrow in fig. [J) is not 
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predicted by the theory, and it cannot be the n = 5 mode which should be at 14.8 ky were its predicted amplitude 
not zero. Obviously, we cannot exclude the possibility that other mechanisms, including Malinkovitch cycles, play 
some role in paleotemperature periodicity. 

Signal emergence times. By the signal "emergence time" we refer to the amount of time required for the nth 
mode amplitude to rise from the modest level of the background noise to a large enough multiple for it to be clearly 
seen in the paleotemperature record. Owing to the large gain, this amplification could occur in a few cycles according 
to the theory. We can determine the emergence time from the data only for modes that turn on after being off for a 
while. Only one of the modes (i.e., n = 2 with a 90 to 100 ky period) clearly shows this behavior, since as seen in fig 
3, it is present only in the first and third time quartiles (solid arrows). It is not easy to tell from fig 1 exactly where 
in the vicinity of 1000 ky before the present this mode begins to appear, and just how sudden is that appearance. 

The shape of the n = 2 cycle, and the abruptness of its emergence can be gleaned by looking at the data for the 
most recent 1000 ky, after making a subtraction of the 41 ky cycle - see fig |U The figure also shows vertical error 
bars on the original data, and a 100 ky cycle drawn to guide the eye. When the data is displayed in this manner the 
100 ky cycle appears much more consistent with a sinusoidal shape than the original data for the last 1000 ky. 

It can easily be seen in fig. 2] that the n = 2 mode appears fairly suddenly (between 800 and 900 ky before the 
present) , following what appears to be a sudden phase reversal or discontinuity. If this discontinuity is real, the theory 
has no obvious explanation for it. However, a short emergence time is consistent with the theory, so this behavior 
offers some support. 

Relative amplitude (and shape) of each mode. 

We have addressed the issue of the lack of a clear signal for the n = 1 mode, and the non-appearance of modes 
with n > 4. Can the theory account for the relative amplitude (and shape) of the three other modes (n = 2,3,4) that 
are clearly seen in the data? It can with the aid of two further assumptions: (1) all modes have the same intrinsic 
strength (as a multiple of the background) when they are active, and (2) only the n = 3 mode is active most of the 
time. 

n = 2 mode (90 ky period). A look at figs.[3]Ja) and (c), i.e., times when both the n = 2 and 3 modes were active, 
shows that the two peaks have comparable strengths. (The position of the n = 2 mode is flagged by solid and dotted 
arrows for those quartiles that it appears to be "on" and "off." ) Thus, the decreased size of this mode in the overall 
data is apparently less a matter of the differences in the intrinsic strength of each mode than what fraction of the 
time this mode is active. The broadening of this peak may also be attributed to its being off perhaps half the time. 
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n — 3 mode (40 ky period). This mode is the dominant (and narrowest) one in the overall record because it is active 
at all times - see figs. Oa) through (d). 

n = 4 mode (22.5 ky period). This mode is not clearly seen in the separate quartile spectra (figs. [3^ a) through (d)). 
However, the fact that it appears to be split in the overall spectrum (into two frequency peaks separated by about 7%) 
has a possible explanation that relates to a single frequency turning on and off in a periodic manner. For example, if 
we perform a Fourier Transform on a sinusoidal signal that turns off and on say every 14 cycles, we would indeed get 
two separate peaks in the frequency domain separated in frequency by 7%. This phenomenon is not due to a defect 
of the Fourier Transform, but rather that two nearby frequencies will produce beats which mimic a single frequency 
turning on and off periodically. It may simply be a coincidence, but this number of cycles on and off for the n = 4 
mode is close to what is seen for the n = 2 mode, which appears to be alternately on and off for approximate intervals 
of 1028 ky (w 11 cycles). 

Finally, how to explain the smallness of the n = 4 peak compared to n = 3? Its splitting into two separate peaks is 
partly responsible, as is the fact that this mode is active only a fraction of the time, unlike the n = 3 mode. A third 
reason for the lower n = 4 amplitude is that the damping of this mode due to blurry boundaries (twice as serious for 
n = 4 than n = 3) operates mainly as the amplification process is occurring while the wave travels through the core, 
which has the effect of preventing the n = 4 mode from reaching a large amplitude. 

Can the oscillations become large enough? 

The dominant (n = 3) mode seen in the Fourier spectrum (fig. 2) is seen to have a signal to background ratio of 
about 3.4. This implies that the background noise at that frequency had to be resonantly amplified by this same 
factor. Given that the complex interactions with the magnetic field that might determine an upper limit to the size of 
the oscillations remains to be worked out, it is unclear if a 3.4-fold growth (at this particular frequency) is achievable. 
This unresolved issue does represent a current shortcoming of the theory. 

VI. PROBLEMS WITH MILANKOVITCH THEORY AND CONCLUSION 

In Milankovitch theory past glaciations are assumed to arise from small quasiperiodic changes in the Earth's 
orbital parameters that give rise to corresponding changes in solar insolation, particularly in the polar regions. A 
brief discussion of five problems with this theory are listed below, and a more detailed description of some of them 
can be found elsewhere. (Karner, 2000) 
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(a) Weak forcing problem: The basic problem with the theory s that observed climate variations are much more 
intense than the insolation changes can explain without postulating some very strong positive feedback mechanism. 

(b) 100 ky problem: The preceding basic problem can be illustrated for the case of one particular parameter 
- the orbital eccentricity. The dominant climate cycle observed during the last million years has a roughly 100 
ky period, which in Milankovitch theory is linked to a 100 ky cycle in the eccentricity. However, the effect of this 
eccentricity variation should be the weakest of all the climate-altering changes, in view of the small change in solar 
insolation it would cause. For example, consider the Earth's orbital eccentricity, e, which has been shown to have 
several periods including one of 100 ky during which e varies in the approximate range: e = 0.03 ± 0.02. (Quinn, 1991) 
The resultant solar irradiance variation found by integrating over one orbit for each of the two extreme e-values is 
about ±0.055%, or ±0.17w/m 2 difference at the top of the Earth's atmosphere. Given that climate models show that 
a one percent change in solar irradance would lead to a 1.8°C average global temperature change, then the change 
resulting from a ±0.055% irradiance change would be a miniscule 0.1°C hardly enough to induce a major climate 
event - even with significant positive feedback. 

(c) 400 ky problem: The variations in the Earth's orbital eccentricity show a 400 ky cycle in addition to the 100 
ky cycle, with the two cycles being of comparable strength. Yet, the record of Earth's climate variations only shows 
clear evidence for the latter. 

(d) Causality problem: Based on a numerical integration of Earth's orbit, a warming climate predates by about 
10,000 years the change in insolation than supposedly had been its cause. 

(e) Transition problem: No explanation is offered for the abrupt switch in climate periodicity from 41 ky to 
100 ky that is found to have occured about a million years ago. Of these five problems with Milankovitch theory, the 
current theory clearly shares only (c). 

In conclusion, We have here suggested a specific mechanism involving diffusion waves in the sun whose amplitude 
should grow very rapidly due to an amplification provided by the link between solar core temperature and luminosity. 
Moreover, the phenomenon of resonant amplification of thermal diffusion waves has been empirically demonstrated, 
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albeit not in the solar context. (Shen, 1995) A number of features of the theory still remain to be resolved, but the 
theory does explain many features of the paleotemperature record, and it appears to be free of most defects of the 
Milankovitch theory. The theory further implicitly suggests the existence of a new category of variable stars having 
extremely long periods - i.e., w 10 4 times longer than stellar periods currently considered to be "very long." For some 
stars with M < M , their thinner radiation zones might make the predicted periods observable. 
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FIG. 1: Proxy record of paleotemperatures back to 5.3 My before present using benthic sediment data. An increase in the 
18 0/ 16 ratio of one part per thousand corresponds to a 5°C decrease in global temperature. 
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FIG. 2: Fourier amplitude spectrum for t = to 4.1 My before present, after subtracting a fifth order polynomial fit to the 
data. The units for the vertical scale are arbitrary. Markers show our theory's predicted periodicities. 
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FIG. 3: Fourier amplitude spectra for (a) to 1023 ky, (b) 1024 to 2047 ky, (c) 2048-3071 ky, (d) 3072 to 4095 ky before 
present. No background subtraction has been made. 
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time in ky before present 

FIG. 4: Residuals for ls O data (including error bars) above a third order polynomial background versus time after making a 
subtraction of the 41 ky (n = 3) contribution, so as to more clearly observe how suddenly the n — 2 mode emerges at 800 ky. 
The data appears to show a phase reversal at that time. 
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TABLE I: Predicted and observed periods (in ky) for modes n = 1 to 4. Uncertainties on the observed values are based on the 
half widths of the peaks in fig. [2] 



mode n 


T n pred 


T„ obs 


1 


360 ± 20 


??? 


2 


90 ±5 


95 ±10 


3 


40 ±2 


41.0 ± 0.4 


4 


22.5 ± 1.5 


(23,7, 22.4) 



